rm(list = ls())

library(parallel)
library(rdrobust)
library(pbapply)
library(tidyverse)

source('code/function_tidy_rd.R')

## Load Civey data with weights 

dat <- readRDS('data/rd_donut_rep.rds')

## Loop over donut BW 

rd_results <- mclapply(dat, function(dfs){
  
  ## Store 
  
  bw <- dfs$donut_bw
  
  ## Loop over data sets (MIP/voting)
  
  pblapply(c(1,2), function(index){
    
    dat_rd <- dfs[[index]]
    
    ## Drop covariate NAs for weighted analysis 
    
    dat_rd_weighted <- dat_rd %>%
      filter(!is.na(weight))
    
    ## Manual bw for weighted RDiT depends on outcome 
    
    bw_manual <- ifelse(index == 1, 18.5, 12.5)
    
    tryCatch({
      
      ##
      ## Estimate RDiT
      ##
      
      ## No weights 
      
      res_rd_unweighted <- rdrobust::rdrobust(y = dat_rd$outcome * 100, 
                                              x = dat_rd$period,
                                              c = 0,
                                              p = 1) %>% 
        tidy_rd(se_nr = 3) %>%
        mutate(bw_exclude = bw,
               outcome = unique(dat_rd$outcome_string),
               w = 'weights: no')
      
      ## With weights (here we use manual bandwidth)
      
      res_rd_weighted <- rdrobust::rdrobust(y = dat_rd_weighted$outcome * 100, 
                                            x = dat_rd_weighted$period,
                                            c = 0,
                                            p = 1,
                                            h = bw_manual, 
                                            weights = dat_rd_weighted$weight) %>% 
        tidy_rd(se_nr = 3) %>%
        mutate(bw_exclude = bw,
               outcome = unique(dat_rd_weighted$outcome_string),
               w = 'weights: yes')
      
      ## Combine and return 
      
      res_rd <- bind_rows(res_rd_unweighted, res_rd_weighted)
      
      
      return(res_rd)
      
    }, 
    
    error = function(e){cat("ERROR :",conditionMessage(e), "\n")})
    
  }) %>%
    reduce(bind_rows) 
}) %>%
  reduce(bind_rows)  %>% 
  mutate(outcome = ifelse(str_detect(outcome, 'Climate'),
                          'Climate change salience',
                          'Green Party support'))

## Plot results 

p1 <- ggplot(rd_results %>%
               filter(w == 'weights: no'), 
             aes(x = bw_exclude, y = estimate, 
                 ymin = conf.low, ymax = conf.high)) + 
  geom_errorbar(width = 0) + 
  geom_point(shape = 21, fill = 'white') + 
  theme_bw() + 
  geom_smooth(se = F, col = 'grey50', method = 'lm',
              formula = y ~ poly(x, 2), size = .5) + 
  facet_wrap(~ outcome, scales = 'free') + 
  geom_hline(yintercept = 0, linetype = 'dotted') +
  scale_x_continuous(breaks = seq(0, 28, 7)) + 
  labs(x = 'Number of days after the flood',
       y = 'Effect estimate (p.p.)') 
p1

## with weights 

p2 <- ggplot(rd_results %>%
               filter(w == 'weights: yes'), 
             aes(x = bw_exclude, y = estimate, 
                 ymin = conf.low, ymax = conf.high)) + 
  geom_errorbar(width = 0) + 
  geom_point(shape = 21, fill = 'white') + 
  theme_bw() + 
  geom_smooth(se = F, col = 'grey50', method = 'lm',
              formula = y ~ poly(x, 2), size = .5) + 
  facet_wrap(~ outcome, scales = 'free') + 
  geom_hline(yintercept = 0, linetype = 'dotted') +
  scale_x_continuous(breaks = seq(0, 28, 7)) + 
  labs(x = 'Number of days after the flood',
       y = 'Effect estimate (p.p.)') 

p2





